My idea for this figure is sort of 1-column wide (~3.5"), with a set of spectral panels going down vertically. The top panel shows a segment of the full resolution model. The next one down can show the individual kernels described above, as well as their composite, $\varphi_v$. Then the next panel down can show the convolved (but still full-resolution) spectrum. Finally, the bottom panel will show the spectrum re-sampled onto the observed pixels. It need not be a large figure (maybe 5" long?).
In [1]:
import matplotlib.pyplot as plt
import numpy as np
from StellarSpectra.grid_tools import HDF5Interface
import StellarSpectra.constants as C
from StellarSpectra.utils import saveall
In [2]:
myHDF5Interface = HDF5Interface("../../libraries/PHOENIX_F.hdf5")
Load the Raw spectrum
In [3]:
wl_raw = myHDF5Interface.wl
fl_raw =myHDF5Interface.load_flux({"temp":6000,"logg":4.5, "Z":0.0, "alpha":0.0})
ind = (wl_raw > 5150) & (wl_raw < 5200)
wl_raw = wl_raw[ind]
fl_raw = fl_raw[ind]
In [4]:
plt.plot(wl_raw, fl_raw)
plt.show()
Plot the kernels in wavelength (velocity) space
In [5]:
@np.vectorize
def gauss_kernel(v, mu=0, FWHM=6.8):
'''
Return the LSF function as a function of velocity
'''
#Convert FWHM to sigma
sigma = FWHM/2.35
return 1./(np.sqrt(2. * np.pi) * sigma) * np.exp(-0.5 * (mu - v)**2/sigma**2)
def vsini_kernel(v, vsini):
'''
Return the vsini function as a function of velocity
'''
epsilon = 0.6
c1 = 2. * (1 - epsilon) / (np.pi * vsini * (1 - epsilon / 3.))
c2 = epsilon / (2. * vsini * (1 - epsilon / 3.))
if np.abs(v) <= vsini:
return c1 * np.sqrt(1. - (v/vsini)**2) + c2 * (1. - (v/vsini)**2) ** 2
else:
return 0
In [6]:
vs = np.linspace(-15, 15, num=200)
k0 = gauss_kernel(vs)
k1 = np.array([vsini_kernel(v, 5.) for v in vs])
#Convolve these two series
k2 = np.convolve(k0, k1, mode="same")
In [20]:
fig = plt.figure(figsize=(3.5, 2))
ax = fig.add_subplot(111)
l0, = ax.plot(vs, k0/np.max(k0), ":", color="0.4", lw=2, label=r"$\mathcal{F}_v^{\rm inst}$")
l1, = ax.plot(vs, k1/np.max(k1), "--", color="0.4", lw=2, label=r"$\mathcal{F}_v^{\rm rot}$")
l2, = ax.plot(vs, 10*np.ones_like(vs), color="0.4", lw=2)
ax.arrow(7., 0, 0, 1, length_includes_head=True, head_width=1.2, head_length=0.08, color="0.4")
l3, = ax.plot(vs + 7., k2/np.max(k2), "k-", lw=2, label=r"$\varphi_v$")
ax.set_xlabel(r"$\Delta v$ [$\textrm{km s}^{-1}$]")
ax.set_ylabel("amplitude")
lg = ax.legend([l0, l1, l2, l3], [r"$\mathcal{F}_v^{\rm inst}$", r"$\mathcal{F}_v^{\rm rot}$",
r"$\mathcal{F}_v^{\rm dop}$", r"$\varphi_v$"], loc="upper left", prop={'size':9})
for txt in lg.get_texts():
txt.set_color("0.2")
fig.subplots_adjust(bottom=0.22, left=0.15, right=0.85)
#fig.subplots_adjust(hspace=0, top=0.93, bottom=0.13, left=0.17)
ax.set_ylim(0,1.1)
ax.set_xlim(-13, 17)
saveall(fig, "../../plots/kernels")
plt.show()
In [11]:
from scipy.interpolate import InterpolatedUnivariateSpline
from StellarSpectra.spectrum import create_log_lam_grid, calculate_min_v
import StellarSpectra.constants as C
from numpy.fft import rfftfreq, rfft, irfft
from scipy.special import j1
In [12]:
FWHM = 6.8
sigma = FWHM / 2.35 # in km/s
vsini = 5
wl_min, wl_max = np.min(wl_raw), np.max(wl_raw)
wl_dict = create_log_lam_grid(wl_min + 4., wl_max - 4., min_vc=0.08/C.c_kms)
wl_FFT = wl_dict["wl"]
min_v = calculate_min_v(wl_dict)
print("min_v is", min_v)
ss = rfftfreq(len(wl_FFT), d=min_v)
ss[0] = 0.01 #junk so we don't get a divide by zero error
interp = InterpolatedUnivariateSpline(wl_raw, fl_raw, k=5)
fl_FFT = interp(wl_FFT)
FF = rfft(fl_FFT)
#Instrumentally broaden the spectrum by multiplying with a Gaussian in Fourier space
taper = np.exp(-2 * (np.pi ** 2) * (sigma ** 2) * (ss ** 2))
#Determine the stellar broadening kernel
ub = 2. * np.pi * vsini * ss
sb = j1(ub) / ub - 3 * np.cos(ub) / (2 * ub ** 2) + 3. * np.sin(ub) / (2 * ub ** 3)
#set zeroth frequency to 1 separately (DC term)
sb[0] = 1.
taper[0] = 1.
FF_tap = FF * sb * taper
fl_blend = irfft(FF_tap, len(wl_FFT))
In [21]:
plt.plot(wl_FFT, fl_blend)
plt.show()
Load TRES data and sample to the same pixels
In [13]:
from StellarSpectra.spectrum import DataSpectrum
mySpec = DataSpectrum.open("../../data/WASP14/WASP14-2009-06-14.hdf5", orders=[22])
wl_TRES = mySpec.wls[0]
ind = (wl_TRES > 5155) & (wl_TRES < 5195)
wl_TRES = wl_TRES[ind]
interp = InterpolatedUnivariateSpline(wl_FFT, fl_blend, k=5)
fl_TRES = interp(wl_TRES)
In [26]:
plt.plot(wl_TRES, fl_TRES, ls="steps-mid")
plt.xlim(5164, 5170)
plt.show()
In [27]:
import matplotlib
from matplotlib.ticker import FormatStrFormatter as FSF
fig, ax = plt.subplots(nrows=2, figsize=(3.5, 3.3), sharex=True, sharey=True)
wl0, wl1 = 5164, 5170
ax[0].plot(wl_raw, fl_raw*1e-7, "r")
ax[0].annotate("raw", (wl1 - 0.2, 1.03), size=9, color="0.2", ha="right")
#ax[1].plot(wl_FFT, fl_blend, "r")
#Doppler shift TRES pixels
wl_TRES_shift = wl_TRES + wl_TRES * (7./C.c_kms)
ax[1].plot(wl_TRES_shift, fl_TRES * 1e-7, "r", ls="steps-mid")
ax[1].annotate("convolved and resampled", (wl1-0.2, 1.03), size=9, color="0.2", ha="right")
ax[1].set_xlim(wl0, wl1)
ax[1].set_ylim(0.0, 1.16)
ax[1].xaxis.set_major_formatter(FSF("%.0f"))
ax[1].set_xlabel(r"$\lambda$ [\AA]")
fig.text(0.01, 0.75, r"$f_\lambda\, \times 10^{7}\,[\textrm{erg}\;\textrm{cm}^{-2}\;\textrm{s}^{-1}\;\textrm{\AA}^{-1}]$",
rotation="vertical")
#for a in ax:
# a.yaxis.set_ticklabels([])
fig.subplots_adjust(hspace=0, top=0.94, bottom=0.13, left=0.15, right=0.85)
saveall(fig, "../../plots/high2low")
plt.show()
In [ ]: